home *** CD-ROM | disk | FTP | other *** search
/ Shareware Grab Bag / Shareware Grab Bag.iso / 098 / fft12.arc / FFT12.PAS < prev    next >
Pascal/Delphi Source File  |  1985-04-26  |  8KB  |  283 lines
  1. program fft;
  2. {$i typedef.sys}
  3. {$i graphix.sys}
  4. {$i kernel.sys}
  5. {$i axis.hgh}
  6. {$i polygon.hgh}
  7.  
  8. type
  9.    cmpl = record
  10.               re,im :real
  11.           end;
  12.  
  13.    bigstr   =  string[12];
  14.    range    =  array[1..256] of cmpl;
  15. var
  16.    filename,nstr,fileout                :bigstr;
  17.    filnam,fiout                         :string[8];
  18.    xform                                :string[25];
  19.    inverse,output                       :char;
  20.    npoints,invt,d                       :integer;
  21.    rd                                   :range;
  22.  
  23. procedure FindWorld(i:integer;
  24.                     A:PlotArray;
  25.                     NPoints:integer);
  26.  
  27.   var XMax,YMax,XMin,YMin:real;
  28.       j:integer;
  29.  
  30.   begin
  31.     NPoints:=abs(NPoints);
  32.     if NPoints>=2 then
  33.       if i in [1..MaxWorldsGlb] then
  34.        begin
  35.         XMax:=A[1,1];
  36.         YMax:=A[1,2];
  37.         XMin:=XMax;
  38.         YMin:=YMax;
  39.         for j:=2 to NPoints do
  40.          begin
  41.           if A[j,1]>XMax then XMax:=A[j,1]
  42.           else if A[j,1]<XMin then XMin:=A[j,1];
  43.           if A[j,2]>YMax then YMax:=A[j,2]
  44.           else if A[j,2]<YMin then YMin:=A[j,2];
  45.          end;
  46.         if (ymax=ymin) and (ymax=0) then begin
  47.            ymax:=1;
  48.            ymin:=0
  49.         end
  50.         else ymin:=0;
  51.         DefineWorld(i,XMin,YMin,XMax,YMax);
  52.         SelectWorld(i);
  53.        end
  54.       else error(7,2)
  55.     else error(7,4);
  56.   end;
  57.  
  58. procedure readdata(filename:bigstr;var d :range;var a:integer);
  59.  
  60. {range and bigstr are set up in the TYPE identifier
  61. range  = array[1..?] of real
  62. bigstr = string[12]
  63. a returns the transform length or file size}
  64.  
  65. type
  66.     data  =
  67.       record
  68.         w,z :real;
  69.       end;
  70. var
  71.    index    :integer;
  72.    xy       :array[1..256] of data;
  73.    datafile :file of data;
  74. begin
  75.      assign(datafile,filename);
  76.      reset(datafile);
  77.      a:=filesize(datafile);
  78.      for index:=1 to a do begin
  79.          read(datafile,xy[index]);
  80.          with d[index] do begin
  81.               re:=xy[index].w;
  82.               im:=xy[index].z;
  83.          end;
  84.      end;
  85.      close(datafile);
  86. end;
  87.  
  88. procedure writedata(filename:bigstr;var d:range;var a:integer);
  89. type
  90.     data  =
  91.       record
  92.         w,z :real;
  93.       end;
  94. var
  95.    index    :integer;
  96.    xy       :array[1..256] of data;
  97.    datafile :file of data;
  98. begin
  99.      assign(datafile,filename);
  100.      rewrite(datafile);
  101.      for index:=1 to a do begin
  102.          with xy[index] do begin
  103.               w:=d[index].re;
  104.               z:=d[index].im;
  105.          end;
  106.      write(datafile,xy[index]);
  107.      end;
  108.      close(datafile);
  109. end;
  110.  
  111. PROCEDURE PLOT(inv,NOPTS :INTEGER; XK :RANGE);
  112.  
  113.       VAR
  114.           HRZ:INTEGER;
  115.           b,a:PlotArray;
  116.           key:char;
  117.           titleRE,titleIM:string[40];
  118.  
  119.    BEGIN
  120.      if inv=1 then begin
  121.         titlere:='Direct Transform --- Amplitude Response';
  122.         titleim:='Direct Transform --- Phase Response';
  123.         FOR HRZ:=1 TO NOPTS DO
  124.           BEGIN
  125.             a[HRZ,1]:=hrz-1;
  126.             b[HRZ,1]:=hrz-1;
  127.             a[HRZ,2]:=sqrt(sqr(XK[HRZ].RE)+sqr(xk[hrz].im));
  128.             b[HRZ,2]:=arctan(XK[HRZ].IM/xk[hrz].re);
  129.           END;
  130.      end
  131.      else begin
  132.         titlere:='Inverse Transform --- Real Part';
  133.         titleim:='Inverse Transform --- Imaginary Part';
  134.         FOR HRZ:=1 TO NOPTS DO
  135.           BEGIN
  136.             a[HRZ,1]:=hrz-1;
  137.             b[HRZ,1]:=hrz-1;
  138.             a[HRZ,2]:=xk[hrz].re;
  139.             b[HRZ,2]:=xk[hrz].im;
  140.           END;
  141.      end;
  142.  
  143.      initgraphic;
  144.      CLEARSCREEN;
  145.      definewindow(1,0,0,xmaxglb,ymaxglb-20);
  146.      definewindow(2,0,0,xmaxglb,ymaxglb-20);
  147.      defineworld(2,b[1,1],-(pi/2),b[nopts,1],(pi/2));
  148.  
  149.      DEFINEHEADER(1,titlere);
  150.      defineheader(2,titleim);
  151.      selectwindow(1);
  152.      setheaderon;
  153.  
  154.      findworld(1,a,nopts);     {automatically executes defineworld}
  155.  
  156.      selectwindow(1);
  157.      drawborder;
  158.  
  159.      drawaxis(7,5,0,0,0,0,0,0,false);
  160.      drawpolygon(a,1,nopts,0,2,0);
  161.      GOTOXY(1,25);
  162.      WRITELN('Press any key to continue...');
  163.      repeat until keypressed;
  164.  
  165.      clearscreen;
  166.      selectworld(2);
  167.      selectwindow(2);
  168.      drawborder;
  169.  
  170.      drawaxis(7,5,0,0,0,0,0,0,false);
  171.      drawpolygon(b,1,nopts,0,2,0);
  172.      repeat
  173.         GOTOXY(1,25);
  174.         WRITE('Press C to continue...');
  175.         read(kbd,key);
  176.      until key in ['c','C'];
  177.  
  178.      leavegraphic;
  179. end;
  180.  
  181. procedure fft(inv,n:integer;var ft:range);
  182.  
  183. {range is set up in the TYPE identifier and is = array[1..?] of real
  184. inv  = 1 if direct transform and =-1 for inverse
  185. n    = transform length power of 2}
  186.  
  187. var
  188.    le,ip,m,i,l,nm1,j,t,le1              :integer;
  189.    ur,ui,tr,ti,tmpr,tmpi,nv2,k          :real;
  190.  
  191. begin
  192.      m:=round(ln(n)/ln(2));
  193.      if abs(m-int(m))=0 then begin
  194.         nv2:=n/2;
  195.         nm1:=n-1;
  196.         j:=1;
  197.         for i:=1 to nm1 do begin
  198.             if i<j then begin
  199.                tmpr:=ft[j].re;tmpi:=ft[j].im;
  200.                ft[j].re:=ft[i].re;ft[j].im:=ft[i].im;
  201.                ft[i].re:=tmpr;ft[i].im:=tmpi;
  202.             end;
  203.             k:=nv2;
  204.             while k<j do begin
  205.                   j:=round(j-k);
  206.                   k:=k/2;
  207.             end;
  208.             j:=round(j+k);
  209.         end;
  210.         for l:=1 to m do begin
  211.             le:=1;
  212.             for t:=1 to l do
  213.                 le:=le*2;
  214.             le1:=round(le/2);
  215.             ur:=1;ui:=0;
  216.             for j:=1 to le1 do begin
  217.                 i:=j;
  218.                 while i<=n do begin
  219.                       ip:=i+le1;
  220.                       tr:=ft[ip].re*ur-ft[ip].im*ui;ti:=ft[ip].im*ur+ft[ip].re*ui;
  221.                       ft[ip].re:=ft[i].re-tr;ft[ip].im:=ft[i].im-ti;
  222.                       ft[i].re:=ft[i].re+tr;ft[i].im:=ft[i].im+ti;
  223.                       i:=i+le;
  224.                 end;
  225.                 ur:=cos(pi*j/le1);ui:=-inv*sin(pi*j/le1);
  226.             end;
  227.         end;
  228.         if inv=-1 then begin
  229.            for i:=1 to n do begin
  230.                ft[i].re:=ft[i].re/n;
  231.                ft[i].im:=ft[i].im/n;
  232.            end;
  233.         end;
  234.      end;
  235. end;
  236.  
  237. begin
  238.      clrscr;
  239.      gotoxy(29,1);writeln('Fast Fourier Transform');
  240.      gotoxy(39,3);writeln('by');
  241.      gotoxy(31,4);writeln('Michael F. Griffin');
  242.      writeln;writeln;writeln;
  243.      gotoxy(19,6);writeln('*----------------------------------------*');
  244.      gotoxy(8,8);writeln('- All file names are assumed to have .DAT as the extension.');
  245.      gotoxy(8,9);writeln('- Output goes to either a file or a printer.');
  246.      writeln;writeln;
  247.      write('Enter data file name w/o extension : ');
  248.      readln(filnam);
  249.      filename:=filnam+'.dat';
  250.      writeln('Input File name ==> : ',filename);
  251.      write('Do you want a direct transform <Y/N>?');
  252.      read(kbd,inverse);writeln;
  253.      if (inverse='Y') or (inverse='y') then begin
  254.         xform:='Direct Fourier Transform';
  255.         invt:=1
  256.      end
  257.      else begin
  258.         xform:='Inverse Fourier Transform';
  259.         invt:=-1
  260.      end;
  261.      write('Do you want the output to a file <Y/N>?');
  262.      read(kbd,output);writeln;
  263.      if (output='Y') or (output='y') then begin
  264.         write('Enter output data file name w/o extension : ');
  265.         readln(fiout);
  266.         fileout:=fiout+'.dat';
  267.         writeln('Output File name ==> : ',fileout);
  268.      end;
  269.      readdata(filename,rd,npoints);
  270.      fft(invt,npoints,rd);
  271.      plot(invt,npoints,rd);
  272.      if (output='Y') or (output='y') then
  273.         writedata(fileout,rd,npoints)
  274.      else begin
  275.         writeln(lst,xform);writeln(lst);
  276.         str(npoints,nstr);
  277.         for d:=1 to npoints do
  278.             writeln(lst,'(',d:length(nstr),')',rd[d].re,rd[d].im)
  279.      end;
  280. end.
  281.  
  282.  
  283.